Main effects are added to the model with a +
DV ~ IV1 + IV2 terms: IV1, IV2
Development and maintenance is independent of “elite” and “esoteric” user groups.
- unlike pretty much all other neuroimaging software tools - R plays nicely with Shell scripts, Python, C++, etc.
FSL and ANTs require Linux (and so does R if you want to use these tools).
When dealing with neuroimaging data in R, we encounter a few important issues with BIG Data
.Machine$integer.max
## [1] 2147483647
# It is good practice to clear your workspace when starting a new project or workflow
rm(list=ls()) # clear all objects from the R environment
gc()# force R to release RAM to the OS
## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 487492 26.1 940480 50.3 940480 50.3
## Vcells 920922 7.1 1650153 12.6 1112749 8.5
Adding packages to your R installation is easy!
install.packages("ggplot2")
pkgs <- c("devtools", "lme4", "lmerTest", "fslr")
install.packages(pkgs)
# Install ANTsR
devtools::install_github("stnava/cmaker")
devtools::install_github("stnava/ITKR")
devtools::install_github("stnava/ANTsR")
# Install TKoscik's Tools
devtools::install_github("TKoscik/tk.nifti") # name change pending / CRAN
devtools::install_github("TKoscik/tk.r.misc") # name change pending / CRAN
devtools::install_github("TKoscik/cluster.nii") # coming soon
devtools::install_github("TKoscik/mRi") # TBD
source("/path/to/your/function.R")
In R your data can take a few different forms
Vector: 1-dimensional object of the same datatype
Matrix: 2-dimensional object of the same datatype
Array: N-dimensional object of the same datatype
Dataframe: 2-dimensional object that allows mixed datatypes
List: 1-dimensional object that contains other objects
R has a suite of utilities to bring your data in
Square brackets, [ ] are used to index your objects
a <- sample(1:10, 5, replace = TRUE)
print(a)
## [1] 8 7 1 2 6
a[3]
## [1] 1
a[c(1,3,5)]
## [1] 8 1 6
df <- data.frame(a = 1:5, b = rnorm(5))
df[2, ]
## a b
## 2 2 0.6388074
df[1,"b"]
## [1] 1.065538
df$a[1:3]
## [1] 1 2 3
temp.ls <- list(a=a, df=df)
temp.ls[[1]]
## [1] 8 7 1 2 6
temp.ls$a
## [1] 8 7 1 2 6
Operators are the symbols that indicate the operation to be implemented
| Operator | Operation |
|---|---|
| <- | leftward assignment to current environment |
| = | leftward assignment (within current context only) |
| -> | rightward assignment |
| <<- or ->> | parent environment assignment |
<- and = are almost interchangeable, however using <- when assigning into a variable the element is sent to the environment as a standalone object as well
| Operator | Operation |
|---|---|
| + | addition |
| - | subtraction |
| * | multiplication |
| / | division |
| ^ or ** | exponentiation |
| %% | modulus |
| %/% | integer division |
| Operator | Operation |
|---|---|
| == | exactly equal to |
| < | less than |
| > | greater than |
| <= | less than or equal to |
| >= | greater than or equal to |
| != | not equal to |
| ! | not |
|| |
OR |
| |
Element-wise OR |
| && | AND |
| & | Element-wise AND |
To do things in R you call functions
* use parentheses ( ) to call a function
* assign output of a function to a variable with <-
* functions have arguments to control what they do
*** Example HERE
One of the most useful concepts in R that is the statistical formula notation that most functions use
~DV ~ IV1 terms: IV1
+DV ~ IV1 + IV2 terms: IV1, IV2
:DV ~ IV1 + IV2 + IV1:IV2 terms: IV1, IV2, IV1:IV2
*DV ~ IV1 * IV2 terms: IV1, IV2, IV1:IV2
: or ’*’DV ~ IV1 * IV2 * IV3 terms: IV1, IV2, IV1:IV2, IV1:IV3, IV2:IV3, IV1:IV2:IV3
DV ~ (IV1 + IV2) * IV3 terms: IV1, IV2, IV1:IV3, IV2:IV3
rm(list=ls())
gc()
# Load ANTs library
library("ANTsR")
Perform bias correction before brain extraction to improve extraction performance. We will redo this using our brain mask to initilize it to get better results after
data.dir <- paste0(getwd(), "/data/debias") # Directory for raw structural data
save.dir <- paste0(getwd(), "/data/debias") # Directory to save processed data
# load raw anatomical image
img.raw <- antsImageRead(filename = paste0(data.dir, "/T1.nii.gz"))
# preform N4 bias correction
img.n4 <- abpN4(img = img.raw,
intensityTruncation = c(0.025, 0.975, 256),
mask = NA,
usen3 = FALSE)
# save N4 corrected image
antsImageWrite(image = img.n4,
filename = paste0(save.dir, "/T1_n4.nii.gz"))
template.dir <- paste0(getwd(), "/data/templates") # Directory for registration and brain extraction templates
# Load brain extraction template
tem <- antsImageRead(paste0(template.dir, "/T_template0.nii.gz"))
# Load registration mask
temmask <- antsImageRead(paste0(template.dir, "/T_template0_BrainCerebellumProbabilityMask.nii.gz"))
# perform brain extraction
img.bex <- abpBrainExtraction(img = img.n4,
tem = tem,
temmask = temmask)
# Save extracted brain
antsImageWrite(image = img.bex[[1]],
filename = paste0(save.dir, "/T1_bex.nii.gz"))
antsImageWrite(image = img.bex[[2]],
filename = paste0(save.dir, "/T1_bmask.nii.gz"))
# preform N4 bias correction
img.n4.bex <- abpN4(img = img.bex[[1]],
intensityTruncation = c(0.025, 0.975, 256),
mask = NA,
usen3 = FALSE)
# save N4 corrected image
antsImageWrite(image = img.n4,
filename = paste0(save.dir, "/T1_n4_bex.nii.gz"))
# Load template image
fixed.img <- antsImageRead(paste0(template.dir, "/MNI152_T1_2mm_brain.nii.gz"))
# Load brain extracted T1
moving.img <- img.n4.bex
# perform rigid, linear, and diffeomorphic (nonlinear) warping
tform <- antsRegistration(fixed = fixed.img,
moving = moving.img,
typeofTransform = "SyN",
initialTransform = NA,
outprefix = paste0(save.dir, "/t1_tform_"))
# apply warping
img.reg <- antsApplyTransforms(fixed = fixed.img,
moving = moving.img,
transformlist = tform$fwdtransforms,
interpolator = "linear",
imagetype = 2)
# save registered image
antsImageWrite(image = img.reg, filename = paste0(save.dir, "/T1_reg.nii.gz"))
func.dir <- paste0(getwd(), "/data/func") # specify directory for raw functional data
func.img <- antsImageRead(paste0(func.dir, "BOLD_run1.nii.gz")) # specify functional data
func.mean <- getAverageOfTimeSeries(func.img)
mocor <- antsMotionCalculation(img = func.img,
fixed = func.mean,
moreaccurate = TRUE,
framewise = TRUE)
mocor.params <- cbind(mocor$moco_params[ ,3:ncol(mocor$moco_params)], mocor$fd)
antsImageWrite(image = mocor$moco_img,
filename = paste0(save.dir, "/BOLD_run1_mocor.nii.gz"))
write.table(x = mocor.params,
file = paste0(save.dir, "/BOLD_run1_mocor.csv"),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = ",")
Register mean BOLD image to structural image
# load fixed image, and resample to 2mm isotropic space
fixed.img <- resampleImage(antsImageRead(paste0(save.dir, "/T1_n4.nii.gz")), rep(2,3))
moving.img <- getAverageOfTimeSeries(paste0(save.dir, "/BOLD_run1_mocor.nii.gz"))
bold.tform <- antsRegistration(
fixed = fixed.img,
moving = moving.img,
typeofTransform = "SyNBold",
initialTransform = NA,
outprefix = paste0(save.dir, "bold_tform_"))
bold.to.t1 <- antsApplyTransforms(
fixed = fixed.img,
moving = moving.img,
transformlist = bold.tform$fwdtransforms,
interpolator = "linear",
imagetype = 3)
bold.to.template <- antsApplyTransforms(
fixed = antsImageRead(paste0(template.dir, "/MNI152_T1_2mm_brain.nii.gz")),
moving = bold.to.t1,
transformlist = tform$fwdtransforms,
interpolator = "linear",
imagetype = 3)
antsImageWrite(image=bold.to.template, filename=paste0(save.dir, "BOLD_run1_norm.nii.gz"))
img.bold <- antsImageRead(paste0(save.dir, "BOLD_run1_norm.nii.gz"))
img.bold.mean <- getAverageOfTimeSeries(img.bold)
img.bold.mask <- getMask(img.bold.mean)
bold.mx <- timeseries2matrix(img.bold, img.bold.mask)
dvars.pre <- computeDVARS(bold.mx)
gc()
compcor.vars <- compcor(fmri = img.bold,
mask = img.bold.mask,
ncompcor = 6,
variance_extreme = 0.975,
randomSamples = 1,
returnv = FALSE,
returnhighvarmat = FALSE,
returnhighvarmatinds = FALSE,
highvarmatinds = NA,
scale = FALSE)
gc()
bold.mx <- residuals(lm(bold.mx ~ scale(compcor.vars)))
gc()
img.bold <- matrix2timeseries(referenceImage = img.bold,
maskImage = img.bold.mask,
timeSeriesMatrix = bold.mx)
antsImageWrite(image = img.bold,
filename = paste0(save.dir, "/BOLD_run1_compcor.nii.gz"))
write.table(x = compcor.vars,
file = paste0(save.dir, "/BOLD_run1_compcor.csv"),
quote = FALSE, row.names = FALSE, col.names = TRUE, sep = ",")
img.bold <- antsImageRead(paste0(save.dir, "/BOLD_run1_compcor.nii.gz"))
img.bold.mean <- getAverageOfTimeSeries(img.bold)
img.bold.mask <- getMask(img.bold.mean)
bold.mx <- timeseries2matrix(img.bold, img.bold.mask)
bold.mx <- frequencyFilterfMRI(boldmat = bold.mx,
tr = antsGetSpacing(img.bold)[4],
freqLo = opts$tfilter.freqLo,
freqHi = opts$tfilter.freqHi,
opt = opts$tfilter.opt)
img.bold <- matrix2timeseries(referenceImage = img.bold,
maskImage = img.bold.mask,
timeSeriesMatrix = bold.mx)
antsImageWrite(image = img.bold,
filename = paste0(save.dir, "/BOLD_run1_tfilter.nii.gz"))
antsImageWrite(image = img.bold,
filename = paste0(save.dir, "/BOLD_run1_prep.nii.gz"))
The temporal sampling of the BOLD signal (one volume at each TR) is usually different from the onsets of your task trials Thus, we want to get the BOLD signal into the space of our task trials. We do this by “Deconvolving” the BOLD signal, which is estimating the strength of the BOLD signal for each task trial.
precision <- 2
params <- list()
params$precision <- 10^precision
params$kernel.length <- 32
params$a1 <- 0.0083333333
params$a2 <- 1.274527e-13
params$p1 <- 5
params$p2 <- 15
params$t <- seq(from = 0, to = params$kernel.length, by = 1/params$precision)
params$hrf <- exp(-params$t) * (params$a1 * params$t^params$p1 -
params$a2 * params$t^params$p2)
params$hrf <- params$hrf[3:length(params$hrf)]
params$hrf <- params$hrf/max(params$hrf)
plot(params$hrf)
ts <- matrix(0,
nrow=(n.tr * tr + params$kernel.length) * params$precision,
ncol = nrow(tf))
for (i in 1:nrow(tf)) {
ts[round(tf$onset[i] * params$precision + 1),i] <- tf$ev[i]
ts[ ,i] <- convolve(ts[ ,i], rev(params$hrf), type = "open")[1:nrow(ts)]
}
conv.ts <- ts[seq(from = 1, to = (n.tr * tr * params$precision), by = tr * params$precision), ]
uf <- data.frame(sim.fmri = rnorm(n.tr, 0,1)*25,
conv.ts)
uf.melt <- melt(uf[ ,c(1,10,20,30)])
uf.melt <- data.frame(1:149, uf.melt)
ggplot(uf.melt, aes(x=X1.149, y=value, color=variable)) +
theme_bw() +
geom_line(size=1)
deconv.mdl <- lm(sim.fmri ~ ., uf)
tf$stbs <- coef(deconv.mdl)[2:(nrow(tf)+1)]
plotf <- melt(tf[ ,c(7,16)])
plotf <- data.frame(trial = 1:36, plotf)
ggplot(plotf, aes(x=trial, y=value, color=variable)) + theme_bw() +
geom_line(size=1)
rm(list=ls())
gc()
library("mRi")
library("parallel")
data.4d <- "/full/path/to/directory/containing/4d/files/or/a/list/of/file/names"
data.paradigm <- "/full/path/to/paradgm.csv"
data.mask <- "/full/path/to/regoin/of/interest/mask(s).nii"
save.dir <- "/empty/directory/to/save/results"
file.prefix <- "some.name"
num.cores <- detectCores()-1
my.command <- "my.function(df, coords, img.dims, save.dir, prefix)"
my.function <- function(df, coords, img.dims, save.dir, prefix) {
mdl1 <- lm(fmri ~ ev, df)
mdl1.coef <- as.data.frame(summary(mdl1)$coef)
table.to.nii(mdl1.coef, coords, img.dims, save.dir, prefix, mdl1)
df$choice <- (df$response + 1) / 2
mdl2 <- glm(choice ~ fmri * ev, df)
mdl2.coef <- as.data.frame(summary(mdl2)$coef)
table.to.nii(mdl2.coef, coords, img.dims, save.dir, prefix, mdl2)
}
run.voxel(data.4d, data.pradigm, my.function, my.command, save.dir, file.prefix, num.cores, vx.order = "rand", verbose = TRUE)
rm(list=ls())
gc()
library("mRi")
library("parallel")
data.4d <- "/full/path/to/directory/containing/4d/files/or/a/list/of/file/names"
seed.vxl <- c(40,40,40)
data.paradigm <- data.frame(seed = load.voxel(parse.4d(data.4d)[[1]], coords = seed.vxl))
data.mask <- "/full/path/to/regoin/of/interest/mask(s).nii"
save.dir <- "/empty/directory/to/save/results"
file.prefix <- "some.name"
num.cores <- detectCores()-1
my.command <- "my.function(df, coords, img.dims, save.dir, prefix)"
my.function <- function(df, coords, img.dims, save.dir, prefix) {
temp <- cor.test(df$fmri, df$seed)
cor.result <- data.frame(r = temp$estimate,
p = temp$p.value,
t = temp$statistic,
df = temp$parameter)
table.to.nii(cor.result, coords, img.dims, save.dir, prefix)
}
run.voxel(data.4d, data.pradigm, my.function, my.command, save.dir, file.prefix, num.cores, vx.order = "rand", verbose = TRUE)
rm(list=ls())
gc()
library(mRi)
cluster.valley(tmap = "",
tvol = 2,
pmap = "",
pvol = 2,
mask = "",
alpha = 0.05,
min.size = 10,
save.dir = "",
file.prefix = NULL)
R has some fantastic visualization tools to make figures ready for publication
rm(list=ls())
gc()
library("mRi")
library("tk.r.misc")
library("ggplot2")
library("gridExtra")
library("raster")
overlay.plot <- draw.overlay(anat.nii = "/full/path/to/anatamoical/image/for/background.nii",
over.nii = "/full/path/to/nii/with/values/to/overlay.nii",
over.vol = 2,
over.color = list(good.color("cool")),
mask.nii = "/full/path/to/mask.nii",
mask.vol = 1,
roi.nii = "/full/path/to/roi/map.nii",
roi.val = 1:10,
roi.color = "#ff64ff",
orientation = "sagittal",
save.plot = FALSE,
return.plot = TRUE)
grid.arrange(overlay.plot, newpage = TRUE)